function new_array = subarray(old_array, add_array, x_offset, y_offset);

% SUBARRAY - add array to portion of old array
%
%    NEW = SUBARRAY(OLD, ADD, X, Y) adds the 2D array ADD to the array OLD
%    at the location specified by (X, Y) and stores the result to NEW.  The
%    size of NEW may be bigger than old if the location & size of ADD require
%    a bigger array.

% By:   S.C. Molitor (smolitor@med.unc.edu)
% Date: February 8, 2000

% initialize output

new_array = old_array;

% validate arguments

if (nargin ~= 4)
   msgbox('Invalid number of arguments', 'SUBARRAY Error', 'warn');
   return
elseif (~isnumeric(old_array) | isempty(old_array) | (length(size(old_array)) ~= 2))
   msgbox('OLD must be a 2D numeric array', 'SUBARRAY Error', 'warn');
   return
elseif (~isnumeric(add_array) | isempty(add_array) | (length(size(add_array)) ~= 2))
   msgbox('ADD must be a 2D numeric array', 'SUBARRAY Error', 'warn');
   return
elseif (~isnumeric(x_offset) | (length(x_offset) ~= 1))
   msgbox('X must be a scalar', 'SUBARRAY Error', 'warn');
   return
elseif (~isnumeric(y_offset) | (length(y_offset) ~= 1))
   msgbox('Y must be a scalar', 'SUBARRAY Error', 'warn');
   return
end
      
% obtain dimensions of NEW
% determine where ADD is to be added

xw = size(new_array, 2);
yh = size(new_array, 1);
xc = [1 size(add_array, 2)] + x_offset - 1;
yc = [1 size(add_array, 1)] + y_offset - 1;

% extend the horizontal dimensions of NEW if range exceeds size

if (xc(1) < 1)
   x_shift = 1 - xc(1);
   new_array(:, x_shift + 1 : x_shift + xw) = new_array;
   new_array(:, 1 : x_shift) = 0;
   xc = xc + x_shift;
   xw = size(new_array, 2);
end
if (xc(2) > xw)
   new_array(:, end + 1 : xc(2)) = 0;
   xw = size(new_array, 2);
end

% extend the vertical dimensions of NEW if range exceeds size

if (yc(1) < 1)
   y_shift = 1 - yc(1);
   new_array(y_shift + 1 : y_shift + yh, :) = new_array;
   new_array(1 : y_shift, :) = 0;
   yc = yc + y_shift;
   yh = size(new_array, 1);
end
if (yc(2) > yh)
   new_array(end + 1 : yc(2), :) = 0;
   yh = size(new_array, 1);
end

% add arrays at specified location

new_array(yc(1) : yc(2), xc(1) : xc(2)) = new_array(yc(1) : yc(2), xc(1) : xc(2)) + add_array;
return
